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The present invention generally relates to apparatus for and 
methods of processing seismic data. It particularly relates 
to methods of performing a decomposition of a seismic 
wavef ield into components such, as up- and downgoing 
wavefield constituents, shear (S) and compressional (P) 
waves and/ or other constituents of interest, where the 
wavefield is obtained through a cross-line survey. 

w 

BACKGROUND OF THE INVENTION 
In the field of seismic exploration, the earth interior is 
explored by emitting low- frequency , generally from 0Hz to 
200 Hz, acoustic waves generated by seismic sources. 
Refractions or reflections of the emitted waves by features 
in subsurface are recorded by seismic receivers. The 
receiver recordings are digitized for processing. The 
processing of the digitized seismic data is an evolved 
technology including various sub-processes such as noise 
removal and corrections to determine the location and 
geometry of the features which perturbed the emitted wave to 
cause reflection or refraction. The result of the 
processing is an acoustic map of the earth interior, which 
in turn can be exploited to identify for example hydrocarbon 
reservo i rs or monitor changes in such reservoirs. 
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Seismic surveys are performed on land, in transition zones 
and in a marine environment. In the marine environment, 
surveys include sources and receiver cables (streamers) 
towed in the body of water and ocean bottom surveys in which 
at least one of sources or receivers are located at the 
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seafloor. Seismic sources and/or receivers can also be 
placed into boreholes. 

The known seismic sources include impulse sources, such as 
5 explosives and airguns, and vibratory sources which emit 
waves with a more controllable amplitude and frequency 
spectrum. The existing receivers fall broadly speaking into 
two categories termed "geophones" and "hydrophones", 
respectively. Hydrophones record pressure changes, whereas 

10 geophones are responsive to particle velocity or 

acceleration. Geophones can recorded waves in up to three 
spatial directions and are accordingly referred to as 1C, 2C 
or 3C sensors . A 4C seismic sensor would be a combination of 
a 3C geophone with a hydrophone. Both types of receivers 

15 can be deployed as cables with the cable providing a 

structure for mounting receivers and signal transmission to 
a base station. 

The spatial distribution of source and receiver locations in 
20 a seismic survey is referred to as layout or spread. A 

variety of spreads are known. Among those are spreads where 
receiver lines, a one-dimensional array of receiver 
locations, and source lines, the corresponding array of 
source or shot locations, are laid out at an angle. For the 
2 5 purpose of this invention, such layouts are referred to as 

"cross-line" geometry or acquisition. Such acquisitions have 
been described for example by G.L.O Vermeer, in U 3D 
Symmetric Sampling", 64th Ann. Internat. Mtg: Soc . of Expl. 
Geophys. (1994), 906-909 and later in the United States 
30 patent no. 6,026,058. 
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Seismic energy acquired at a receiver may contain upwardly 
and/or downwardly propagating seismic energy depending on 
the location of the receiver and on the event. For example 
seismic energy when it is incident (travelling upwardly) on 
5 the water-seabed interface, be partly transmitted into the 
water column and partially reflected back into the seabed. 
Thus, a seismic event will consist purely of upwardly 
propagating seismic energy above the seafloor, but will 
contain both upwardly and downwardly propagating seismic 

10 energy below the seafloor. As another example, seismic 

energy when incident on the water-air interface at sea level 
will be reflected back into the water column generating so- 
called "ghost" events. It is therefore often of interest to 
decompose the; seismic data acquired at the receiver into an 

15 up-going constituent and a down-going constituent. 

Various filters that enable decomposition of seismic data 
into up-going and down-going constituents have been 
proposed. For example in w Application of Two-Step 
2 0 Decomposition to Multi-Component Ocean-Bottom Data: Theory 
and Case Study", J . Seism. Expl . Vol. 8, 261-278 (1999), 
K.M. Schalkwijk et al have suggested that the down- going and 
up-going constituents of the pressure just above the 
seafloor may be expressed as: 

25 

[1] 

P~(f, k x , ky) = I P^, k x , k y ) ~ ^ ^ V^f, k x , k^ 

P^f, k*, ky) = l Kf. k x , ky) + 2q(£# ^ # ^ ky) ' 

where P is the pressure acquired at the receiver, P T is the 
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up-going constituent of the pressure above the seafloor, P + 
is the down-going constituent of the pressure above the 
seafloor, f is the frequency, k x , k y are the horizontal 
wavenumbers, v 2 is the vertical particle velocity component 
5 acquired at the receiver, p is the density of the water, and 
q is the vertical slowness in the water layer. 

As can be seen, the expressions in equation [1J require two 
of the components of seismic data recorded at the receiver 

10 to be combined. These expressions are examples of combining 
two components of the acquired seismic data. It may also be 
necessary to combine two or more components of the acquired 
seismic data in order to decompose the acquired seismic data 
into P-wave and S-wave components, or to remove water level 

15 multiple events from the seismic data. 

Further separation methods including free-surface multiple 
removal above the seafloor, wavefield decomposition into up- 
and downgoing constituents or P/S events above and below the 
20 surface, the splitting of particle velocities and traction 
are described in a number of published documents. 

In United States patent no. 6,101,408, the ocean bottom 
wavefield separation described in three dimensions using an 

25 analytical solution. However, for practical applications, 
the filter is reduced to one dimension. A number of 
decomposition equations for various separations are 
developed by Amundsen et al. in the above cited United 
States patent no. 6,101,408 and in: "Multiple attenuation 

30 and P/S splitting of multicomponent OBC data at a 

heterogeneous sea floor". Wave Motion 32 (2000), 67-78- A 
further review of decomposition methods for use in 
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connection with the present invention is presented by Xj- 
Amundsen in: -Elimination of free-surface related multiples 
without need of the source wavelet", Geophysics, Vol. 66, 
No. 1 (Jan-Feb 2001), 327-341. 

Approximated compact spatial filters are further described 
by Osen et al. in: Towards Optimal Spatial Filters for 
Multiple Attenuation and P/S-Splitting of OBC Data", EAGE 
60 th conference, Leipzig, Germany, 8-12 June 1998, 1-29 
Geophysical Division. A short length filter is obtained in 
terms of powers of k x using a series expansion. 

When applying three-dimensional (3D) wavefield decomposition 
methods to data acquired in a cross-line geometry and sorted 
into 1-fold bins of common mid-points (CMPs) distributed 
evenly in a finely spaced "carpet" determined by in-line 
source and receiver spacings as proposed by Vermeer, it was 
noted that the known filter introduce an unacceptable level 
of noise due to sensor variations, statics and other 
perturbations . 

In the light of the above prior art, it is seen as an object 
of the present invention to provide filters applicable to 
cross-line acquisitions or data collected through cross-line 
acquisitions and methods of applying such filters. 

SUMMARY OF THE INVENTION 

According to a first aspect of the invention there is 
provided a method of decomposing a seismic wavefield, 
wherein a 3D wavefield is obtained by a cross-line 
acquisition and filtered applying a decomposition filter 
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having two spatial components or filtering in two spatial 
directions to obtain a decomposed representation of the 
wavef ield. 

A 3D wavefield for the purpose of this invention involve 
obtaining data or time series of measured parameters over an 
area. Hence, such data are acquired as series of ideally 
closely spaced parallel lines. The parameters measured are 
preferable velocity and pressure data. 

In the cross-line acquisition stpt lines and receiver lines 
enclose an angle, which is preferably around 90 degrees. 

The method of this invention can be applied to any of the 
existing decomposition equations that include a filter term 
depending on the vertical wavenumber k 2 . Such decompositions 
preferable include up- / down going decomposition, P/S 
decomposition, elastic decomposition and acoustic 
decomposition . 

The filter of the present invention has two spatial 
components that when represented in an analytical form are 
written as k x and k y , or as spatial derivatives in x and y, 
respectively- When implemented as machine program these 
filters are approximated by finite differences. 

The filter is preferably a cascaded filter of ID spatial 
filters that are applied sequentially. 

The filter is • preferably a compact filter having a finite 
length or support in in-line and cross-line direction. When 
analytically derived, wavefield decomposition filters have 
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infinite extent or support in space and in time. The filter 
operations generally assume stationary medium properties or 
in the case of deghosting a locally flat sea surface (both 
in time and space) . The main advantage of introducing 
5 compactness is to ensure that medium properties (or sea 

surface variations) are constant across the aperture (or two 
dimensional support) of the filter. For the seabed 
wavefield decomposition filters the full analytical 
expression of the filter can be written in the frequency 

10 wavenumber domain for instance. The infinite support in 

time can be maintained since medium properties do not vary 
with time. However, in the spatial directions Taylor 
Gxpemslons of the filter into factors k x/ k x ^2, k x A 3, 
k y/ k y A 2, k y A 3, ... are proposed. When going back to the 

15 spatial domain each factor k x or k y or its powers simply 
correspond to a derivative in the x- and y-directions 
respectively. Spatial derivatives can in turn be 
implemented with compact local support using 2 -point, 3- 
point, 5 -point, or more extended FD approximations . 

2 0 However, there are also other ways of designing compact 
filters without necessarily relating them to spatial 
derivatives. 

Preferably the spatial filter of the present invention is 
25. applied exclusively to the measured pressure wavefield 
P(f,k x ,k y ). The pressure measurement is usually less 
sensitive to mismatch in the response of the various 
receivers used to record the wavefield. 

30 It is furthermore advantageous to use a calibration for 
matching geophone recordings with hydrophone recordings 
prior to the decomposition filtering, particularly in case 
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the filter operates on the particle velocity (v x/ v y ,v z ) 

After the decomposition filter is applied it is possible to 

remove multiples or proceed with other known steps to obtain 

5 an image of sub-surface/ including migration and other 
methods known in the art . 



These and other aspects of the invention will be apparent 
from the following detailed description of non- limitative 
10 examples and drawings, 

BRIEF DESCRIPTION OF THE DRAWINGS 
FIG. 1A illustrates an ocean bottom acquisition of a 3D 

seismic wavefield using an OBC and a source towed 
15 by a seismic vessel; 

FIG. IB illustrates the spread of data point after shooting 

a single shot line of the acquisition of FIG. 1A; 

2 0 FIG. 1C illustrates the acquisition of 3D seismic wavefield 

using streamers towed by a first vessel and a 
source towed by a second seismic vessel; 

FIG. 2 is a diagram illustrating steps in accordance with 
25 an example of the invention; and 

FIG. 3 compares the performance of filter approximations s 

in accordance .with examples of the invention. 



3 0 EXAMPLES 

In FIG.1A there is illustrated an example of a seismic 
survey in cross-line geometry. The survey is a marine 
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seismic survey in a body of water 10 between a seafloor 101 
and a sea surface 102. A receiver cable 11 with a plurality 
of receivers 111 is laid out on the seafloor 101. The 
receivers 111 are preferably 4C sensors, though, as will be 
5 apparent from the following description, the sensors may 
comprise a 1C geophone and a hydrophone only, being thus 
capable of recording at least the vertical component of 
velocity and pressure at seafloor level. 

A seismic vessel 12 tows a marine seismic source 13 close to 
the sea surface 102. The airgun 13 emits at precisely 
determined time intervals an impulse of acoustic energy 
referred to as "shot". A dashed line 132 indicates the path 
of the towed airgun 13. The projection 133 of the dashed 
line 132 onto the seafloor 102 intersects the receiver line 
112 at approximately 90 degrees. Though it is preferable to 
aim for a near-orthogonal orientation of receiver lines to 
shot lines, deviations are inevitable under real survey 
conditions. To facilitate the following description the 
receiver line or in-line direction is denoted as x 
direction, the shot-line or cross-line direction is marked 

« 

as y direction and the vertical direction is taken as the z 
direction . 

25 During a survey, the sources 131 are fired at intervals and 
the receivers 121 "listen" within a frequency and time 
window for acoustic signals such as reflected and/or 
refracted signals that are caused features in path of the 
emitted wavefield. After shooting a line, the vessel 

30 performs a u-turn in order to shoot a subsequent line 

usually with an offset in receiver line or x direction. 
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In the general practice, it is assumed that Green's 
functions which describe the wave propagation between source 
and receiver points are invariant for translation of the 
source and receiver in the cross-line direction. Hence/ an 
5 offset between shooting lines can be regarded as an equal 

■ 

shift of the receiver line. As a result, data points 
obtained from a single source cross-line form a carpet on 
the seafloor, which is illustrated in FIG . IB. In FIG . IB, 
the triangles 122 denote the location of data points. As 
10 the wavefield is recorded in two spatial dimensions and in 
time, the resulting data are referred to as 3D wavefield. 

In FIG. 1C, there is shown a schematic cross-line marine 
survey with two vessels. A first vessel 15 tows five 
15 streamers 151 below the sea surface following a path 152. 
Simultaneously, a second vessel 16 on path 162 tows a 

seismic source 161 below the surface. As above in Fig- 1A, 

> 

the resulting shot and receiver lines are essentially 
orthogonal to each other. 

20 

In the following description and the accompanying FIG. 2, 
steps are described leading to a decomposition of the 3D 
wavefield into up- and downgoing components. 

25 After obtaining 20 the wavefield data set as acquired 

through the seismic receivers, the data are first preferably 
calibrated 21 to compensate for the differences between 
geophone and hydrophone recordings. Any suitable 
calibration may be used including for example the methods 

30 described in the International patent application 

PCT/GB03/04190 . Following those methods, the calibration 
can be done using an in-line shot- line for the P, v z and v z 
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components and using a cross- line shot-line for the v y 
component . 

Acoustic wavefield decomposition 22 is usually carried out 
5 on the pressure component P (involving spatial filtering of 
v z ) . Instead in this example decomposition filters are 
applied to the vertical geophone component v z {involving 
spatial filtering of P) . The advantage of this example is 
that the spatial components of the filter only act on P as 
10 shown in PCT/GB03/ 04190 , 

Accordingly, acoustic wavefield decomposition into up- and 
down-going constituents above the seafloor can be achieved 
by solving the following equation: 



15 



[2] 



vr <f, k x , k^ = - a (f ) v z (f, k x , ky) ± — ^- P (f, k x , ky) 

z x -y 2 y 2pco * 

In equation [2], a(f) denotes the optional frequency- 
20 dependent calibration filter that corrects for imperfections 
in the recording of the geophones, v z " denotes the up-going 
constituent of the vertical component of particle velocity, 
and v z + denotes the down-going constituent of the vertical 
component of particle velocity. The velocity v z is the 
25 recorded or estimated vertical component of particle 

velocity in the frequency f - wavenumber domain, P(f,k X /k y ) 

is the recorded pressure, and p is the density in the 
recording medium. 
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The term k z , which can be expressed as 
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[3] k 2 (f,k X# ky) = ^Ttf/C) 2 - - ky 2 , 

is the absolute value of the vertical wavenumber expressed 
in terms of horizontal wavenumbers in the in-line direction 
5 k x and the cross-line direction k y , and the velocity c of 
the recording medium. It should be noted that the 
decomposition could also be achieved by computing the up- 
going component of the recorded pressure P~ using equations 
[1] , leading expression which include terms of l/k z . Such 
10 terms can be approximated using similar expansions as 
described below. 

In known decomposition methods using any of the above 
equations [1,2], the cross-line or y-directions is mostly 
15 ignored or a radial symmetry is assumed, with the vertical 
wavenumber then being computed using an approximation based 
exclusively on a one-dimensional direction, i.e. the in-line 
wavenumber k x or the radial wavenumber k r . When 3D effects 
of the sub- surface or acquisition geometry are significant, 

20 such approximations are no longer valid. 

-. 

Equations as those described herein can be implemented in 
the common mid-point domain, which is proposed by Vermeer 
(1994) and Thomas (2000) . It is however fruitful to rewrite 

25 or approximate equation [3] into a form constituting a 

cascade (sum or product) of one-dimensional (ID) spatial 
filters acting in the x- or y-directions only. This 
represents a computationally attractive way of filtering the 
data (both in terms of CPU and resorting data between 

30 different domains) • One way to obtain filters of this form 
is to make suitable Taylor expansions of the horizontal 
wavenumbers in the square-root term around zero wavenumbers . 
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This approximation remains valid for data corresponding to 



propagating waves at k x J + ky^ < {2n£/c) z ) . 



The expression for the vertical wavenumber can be rewritten 
and expanded in k x and k y to produce a few different 
alternative expansions that can be implemented using a 
cascade of filters that only act in the cross-line or in- 
line direction one at a time: 



10 [4a] 



k z (f/ k x # — 



A /(27tf/c) 2 - k x 2 



1 - 



V 



+ 



*y 2 

2((27tf/c) 2 - k x 2 ) 
o(k x 6 ,ky 6 ) 



8((27if/c) 2 - k x 2 f 



J 



15 



[4b] 



k 2 (f # k x jr k v ) — 



( 



2n£/c 



1 - 



V 



2(2rcf/c) 2 



k x + ky + 2k x ^ky 2 



\ 



8(27if/c) 4 



+ 



O^ky 6 ) 



J 



[4c] 
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k 2 (f* ICj^r ky) ~ 



• > /(27tf/c) 2 - k x 2 + - N /(27if/c) 2 - ky 2 - 27tf/c 



1 i ^x 2 ^ 

8(2«f/c) 4 
v + o(k x 6 ,k y 6 ) y 



Equations [4a-4c] represent different ways to proceed with 
an implementation of the filter with two spatial components 
5 that only rely on being able to filter the data along two 
perpendicular spatial directions one at a time which is 
exactly what can be achieved using the method described 
above in a cross-line geometry. Note that this does not 
mean that only a *cross" of mid-points are used in filtering 
10 the data. All cross-terms of multiplications of terms with 
different horizontal wavenumbers will result in a "virtual" 
carpet of data being used of dimension of the length of the 
spatial filters in both directions. 

* 

15 After the decomposition 22, multiples could be removed 23 
from the data set. Further processing steps and/or 
filtering steps 24 could be performed on the decomposed data 
set. Using what is commonly referred to as imaging or 
migration 25 the data set can be further processed to yield 

2 0 an image of subterranean layers. These images are used for 
hydrocarbon exploration and reservoir characterization. The 
optional steps 21 and 23 - 25 are indicated in FIG. 2 as 
dashed blocks . 

25 FIG - 3 shows a panel of the exact wavefield decomposition 
filter using equation [3] in the top left, difference 
between equation [3] and the filter approximation [4a] in 
the top right, difference between equation [3] and the 
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filter approximation [4b] in the bottom left, and difference 
between equation [3] and the filter approximation [4c] in 
the bottom right. Note that these plots only assess the 
accuracy of the filter approximations and do not include the 
5 error due to their discretization. In other words, the 

plots do not show inaccuracies related to how the different 
terms in the spatial filter approximations [4a-4c] are 
implemented (e.g., using 3 -point or 5 -point derivative 
approximation) . This would of course introduce a dependence 
10 on frequency as well. However, this is of secondary 
importance as appropriate approximations that are 
sufficiently accurate are straightforward to find. 

From the top right of FIG* 3, it can be seen that equation 

15 [4a] which was obtained by making a Taylor expansion in the 
y-direction only results in the best approximation of the 
three examples for azimuths close to the in-line cable 
direction. Equation [4b] which is used in the difference 
plot at the bottom left of FIG. 3 results in an 

20 approximation which is equally good along all azimuths. An 
advantage with this filter is that it can be fully 
implemented as a compact filter. Equation [4c] which is 
used in the difference plot at the bottom right of FIG. 3 
results in a fully accurate approximation both along the in- 

25 line and cross-line azimuths. The filter can be implemented 
using a compact filter approximation for the cross-term 
only. Exactly which of the alternative implementation [4a- 
4c] that is most attractive may vary depending on different 
combinations of requirements in terms of computational cost 

3 0 (CPU and data access) and accuracy. 
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While the invention has been described in conjunction with 
the exemplary embodiments described above, many equivalent 
modifications and variations will be apparent to those 
skilled in the art when given this disclosure. Accordingly, 
5 the exemplary embodiments of the invention set forth above 
are considered to be illustrative and not limiting. Various 
changes to the described embodiments may be made without 
departing from the spirit and scope of the invention. 

10 The above approximations or similar approximations can be 

used for example with the separations developed by Amundsen 
et al. in the above cited United States patent no. 6,101,408 
or in: "Multiple attenuation and P/S splitting of 
multi component OBC data at a heterogeneous sea floor", Wave 

15 Motion 32 (2000), 67-78. In the latter document, demultiple 
or decomposition equations are found for elastic 
decomposition (particle velocity, traction) or P/S wave 
splitting below the sea floor. 
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